R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/oisst_an.m
First calculate 1982-2011 SST monthly mean. To evaluate SST anomaly:
SST.1 = Monthly SST - SST_30yr_mean
SST.anomaly = SST.1 - Linear trend by monthly series from 1982-202205
Land mask from R/01-2_OISST_land_seaice_mask.Rmd
## The same code from R/01-2_OISST_land_seaice_mask.Rmd(and can check that testing figure)
land_mask <- fread("unzip -p ../data_src/mapfiles/sea_icemask_025d.csv.zip")[, .(lon, lat, x, y, landmask)]
setorder(land_mask, -y, x)
land_mask[,ry:=rowid(x)]
yy <- unique(land_mask[,.(y, ry)]) ## check, NOTE that in OISST.nc file, y is from Northest (89.875) to Southest (-89.875)
landm <- dcast(land_mask, x ~ ry, value.var = "landmask")
xx <- as.numeric(landm[,1]$x)
landm <- as.matrix(landm[,-1])
clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
Recalculate_30yrs = FALSE
if (Recalculate_30yrs) {
plan(multisession)
options(future.globals.maxSize= 1048576000)
stm <- future_lapply(1:12, function(j) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in clim_years) {
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- jstr
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_
if (i == clim_years[1]) {
stmx <- x
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stmx <- c(stmx, x)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
names(stmx) <- rep(jstr, length(names(stmx)))
stmx <- merge(stmx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>%
aggregate(by=paste0(climyrs, " years"), FUN=mean, na.rm=TRUE)
names(stmx)[1] <- jstr
stmx <- stmx %>% select(jstr) %>% adrop
return (stmx)
})
}
if (Recalculate_30yrs) {
# Note that in an earlier version used "GLDASp4 land mask" but current version use "GLDASp5"
# The differnt version of land mask cause difference in the definition of land areas, so that sea-mask
# and then also cause difference in climatology mean result
save(stm, file="../data_src/stats/sst_1982_2011month_climatology_nc.RData")
} else {
load("../data_src/stats/sst_1982_2011month_climatology_nc.RData")
}
stmlibrary(grid)
plotx <- vector("list", length = 12)
pstrx <- '';
for (j in 1:12) {
plotx[[j]] <- gplotx(stm[[j]], monstr[j], returnx = TRUE, minz = -2.5, maxz = 34.5)
pstrx <- paste0(pstrx, "plotx[[", j, "]],")
}
lay1 <- rbind(c(1,1,2,2,3,3),
c(4,4,5,5,6,6),
c(7,7,8,8,9,9),
c(10,10,11,11,12,12))
evplot <- paste0('grid.arrange(', pstrx, ' layout_matrix=lay1)')gx <- eval(parse(text=evplot))To double check if the plots of climatology are right, use existed data that download from NOAA (NetCDF OISST, 0.5d, monthly mean)
Metadata: https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html
https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc (downloaded 202206 and the data is 1981/12 to 2022/05)
if (!require("ncdf4")) install.packages("ncdf4"); library(ncdf4)
sst_mnfile = "../data_src/stats/sst.mnmean.nc"
if (!file.exists(sst_mnfile)) {
tryCatch({
curl::curl_download("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc", destfile = sst_mnfile)
}, error = function(e) paste0(e))
}
nx0 <- nc_open(sst_mnfile)
print(nx0) ## sst[lon,lat,time] (Chunking: [360,180,464])
latn1<- ncvar_get(nx0, "lat") # 89.5 - -89.5
lngn1<- ncvar_get(nx0, "lon") # 0 - 359.5
time<- ncvar_get(nx0, "time") # 1981-12.01 - 2020-07-01, we now check every March, August, and December
date<- time %>% as.Date(origin="1800-01-01 00:00:0.0")
maridx <- seq(4, 464, by = 12) #check date[maridx]
augidx <- seq(9, 464, by = 12)
decidx <- seq(13, 464, by = 12)
get_sstmx <- function (midx, minz=-2.5, maxz=34.5) {
dtx <- as.data.table(ncvar_get(nx0, "sst")[,,midx]) %>%
.[,.(sstm=mean(value, na.rm=TRUE)), by=.(V1,V2)] %>%
.[,`:=`(longitude=fifelse(lngn1[V1]>180, lngn1[V1]-360, lngn1[V1]), latitude=latn1[V2])] %>%
.[,.(longitude, latitude, sstm)]
setnames(dtx, 1:2 , c("x", "y"))
gx <- gplotx(dtx, var="sstm", returnx=TRUE, minz=minz, maxz=maxz, no_coord=TRUE, maxlim=180)
gx <- gx +
#gx <- ggplot() +
# geom_tile(data=dt, aes_string(x="longitude", y="latitude", fill="sstm"), alpha=0.8) +
# scale_fill_viridis(limits=c(minsst, maxsst)) +
geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) #+
# coord_sf() +
# xlim(c(-180, 180)) + ylim(c(-90, 90))
return(gx)
}
gx3 <- get_sstmx(maridx)
gx8 <- get_sstmx(augidx)
gx12 <-get_sstmx(decidx)gt3 <- gplotx(as.data.table(stm[[3]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Mar", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
gt8 <- gplotx(as.data.table(stm[[8]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Aug", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
gt12 <- gplotx(as.data.table(stm[[12]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Dec", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
lay2 <- rbind(c(1,1,2,2),
c(3,3,4,4),
c(5,5,6,6))grid.arrange(gx3, gt3, gx8, gt8, gx12, gt12, layout_matrix=lay2)i <- 2022
j <- 4
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- "anom"
#dim(x[[1]]) #1440. 720
#dim(landm) #1440, 720
#dim(ice) #1440, 720
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_
#xt <- matrix()
#length(xt) <-1440 * 720
#dim(xt) <- c(1440, 720)
## if one of x and stm is NA, then the substraction is NA (so numbers of NA values increases)
x[[1]] <- x[[1]] - stm[[j]][[1]]## Just check if plot 2022-04
gx <- gplotx(x, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
names(y)[1] <- "anom"
gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)
layt <- rbind(c(1,1,2,2))
grid.arrange(gx, gy, layout_matrix=layt)yrng <- seq(1982,2022)
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))
j=4 #just testing, arbitrary select a month
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in yrng) {
if (i==curryr & j>(currmo-1)) break
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
if (i==yrng[1]) {
stdrx <- z
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stdrx <- c(stdrx, z)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
names(stdrx) <- rep(jstr, length(names(stdrx)))
xt <- merge(stdrx)
st_set_dimensions(xt, 3, values = as.POSIXct(datex), names = "time")## stars object with 3 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## X -2.8777 -0.4862444 -0.09030438 -0.2813758 0.05061904 1.245278 95907
## dimension(s):
## from to offset delta refsys point
## x 1 1440 0 0.25 NA NA
## y 1 720 90 -0.25 NA NA
## time 1 41 NA NA POSIXct NA
## values x/y
## x NULL [x]
## y NULL [y]
## time 1982-04-01 08:00:00,...,2022-04-01 08:00:00
seqx <- rbind(c(`377`= as.numeric(xt[[1]][377,41,])), ## arbitrary select a point, check which(!is.na(xt[[1]][,41,]))
c(`1112`= as.numeric(xt[[1]][1112,41,])),
c(`test`= sort(rnorm(41)) + rnorm(41)))
colnames(seqx) <- gsub("-", "_", datex)
rownames(seqx) = c("0377_41", "1222_41", "test")
yhat <- stats::predict(lm(seqx[3,] ~ seq_along(datex)))
seqx[1, 28:30] <- NA_real_ #insert NA
seqx[2, 32:34] <- NA_real_ #insert NA
predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
haty <- apply(seqx, MARGIN=1, FUN=function(y){ #predy use predict will auto exclude NA, so that out length is shorter
if (all(is.na(y))) return(y)
yh <- as.numeric(predy(y))
y[!is.na(y)] <- yh
return(c(as.numeric(y)))
}) %>% as.data.frame()
all.equal(haty$test, as.numeric(yhat)) #TRUE## [1] TRUE
#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)## Loading required package: pracma
##
## Attaching package: 'pracma'
## The following objects are masked from 'package:magrittr':
##
## and, mod, or
ydet <- as.numeric(detrend(seqx[3,], 'linear'))
all.equal(ydet, as.numeric(seqx[3,]) - as.numeric(yhat)) #TRUE## [1] TRUE
#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)
plot_pred <- function(x, y, yhat=NULL) {
if (is.null(yhat)) {
yhat <- stats::predict(lm(y ~ x))
} else {
yh <- stats::predict(lm(y ~ x))
if (!isTRUE(all.equal(as.numeric(yh), yhat))) {
print(paste0("Not equal yhat with lm: ", paste(yh, collapse=",")))
}
}
plot(x=x, y=y)
abline(lm(y ~ x))
points(x=x, y=yhat, col="red", pch=4)
points(x=x, y=c(y-yhat), col="green", pch=13)
}
plot_pred(seq_along(datex), seqx[3,], haty$`test`)save(seqx, file="../data_src/stats/test_stars_predict00.RData")Recalculate_anom <- FALSE
if (Recalculate_anom) {
for (i in yrng) {
mmx <- fifelse(i==curryr, currmo-1L, 12L)
for (j in 1:mmx) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
filet <- paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc")
if (!file.exists(filet)) {
#z<- read_stars(paste0("../data_src/oisst/monthly_anom/", i, monj, "_anom.nc")) #Old vers use anom, not sst, so no subtract 30yr mean
z <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(z)[1] <- "anom"
z[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
z[[1]][which(ice[[1]]==1)] <- NA_real_
z[[1]] <- z[[1]] - stm[[j]][[1]] #Old-version without this substraction while using _anom directly
write_stars(z, filet)
}
}
}
}## Just check if plot 2022-04
#if (Recalculate_anom) {
gx2 <- gplotx(z, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
#y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
#names(y)[1] <- "anom"
#gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)
layt <- rbind(c(1,1,2,2))
grid.arrange(gx2, gy, layout_matrix=layt)#}## It's hard/slow to detrend each x,y along almost 30yr (30x12) year trend. Must read 1440x720x30x12 at once...
#removeLargeObj <- TRUE
#if (removeLargeObj) {
# rm(land_mask, landm, stm)
#}
# 20220727 modified for detrend not distinguish month, all time-span should involved
Reload_allnc <- FALSE
#predyr <- 2021 # previous year that has whole 1-12 month data
#predrng<- seq(1982, predyr)
yrs <- yrng[length(yrng)] - yrng[1] + 1
#yrs <- predrng[length(predrng)] - predrng[1] + 1
if (Reload_allnc) {
for (i in yrng) {
for (j in 1:12) {
if (i==curryr & j>(currmo-1)) break
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
names(z) <- "anom"
if (i==yrng[1] & j==1) {
stdrx <- z
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stdrx <- c(stdrx, z)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
}
print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
# save(stdrx, file="../data_src/stats/monthly_anom_1982_2021_nc.RData") #too big, not save
} #else {
# load("../data_src/stats/monthly_anom_1982_2021_nc.RData")
#}Recalculate_anom <- FALSE
if (Reload_allnc) {
trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time")
tx <- st_get_dimension_values(trd, 3)
predtx <- function(x) {
if (all(is.na(x))) return(list(as.numeric(x)))
x[!is.na(x)] <- as.numeric(stats::predict(lm(x ~ tx)))
return(list(as.numeric(x)))
}
trdx <- st_apply(adrop(trd), c(1,2), predtx)
tk <- array(unlist(trdx[[1]]), dim = c(length(datex),1440,720))
save(tk, file="../data_src/stats/predict_anom_1982_202204_nc.RData") #too big, not save
} else {
if (Recalculate_anom) load("../data_src/stats/predict_anom_1982_202204_nc.RData")
}if (Recalculate_anom) {
for (k in seq_along(yrng)) {
for (j in 1:12) {
if (yrng[k]==curryr & j>(currmo-1)) break
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
nt <- (k-1)*12+j
# for (i in yrng) {
# if (i==curryr & j>(currmo-1)) break
#
# z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
#
# if (i==yrng[1]) {
# stdrx <- z
# datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
# } else {
# stdrx <- c(stdrx, z)
# datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
# }
# }
# names(stdrx) <- rep(jstr, length(names(stdrx)))
# print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
# predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
# trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>%
# aggregate(by=paste0("months"), FUN = function(y) {
# if (all(is.na(y))) return(list(as.numeric(y)))
# y[!is.na(y)] <- as.numeric(predy(y)) #Seems if return value has equal length, the result will be simplified by aggregate if using as.matrix
# return(list(as.numeric(y))) #Old-version bug: not predy(y), should be y. See previous testing chunk
# })
# tk <- array(unlist(trd[[1]]), dim = c(length(datex),1440,720)) #bug fix: NOT each predict have the same length result (if input has NA)
# print(paste0("Predict ok with dim(tk): ", paste0(dim(tk), collapse=",")))
# for (k in seq_along(yrng)) {
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", yrng[k], monj, "_anom.nc"))
#notna1 <- which(!is.na(z[1][[1]]))
#notna2 <- which(!is.na(tk[k,,]))
#notna <- intersect(notna1, notna2)
#z[[1]][notna] <- z[[1]][notna]- tk[1,,][notna]
z[[1]] <- z[[1]]- tk[nt,,] #Any of z[[1]] or tk[1,,] NA will cause NA here, it makes sense.
names(z)[1] <- "anom"
filet <- paste0("../data_src/oisst/monthly_anom_detrend/", yrng[k], monj, "_anom.nc")
write_stars(z, filet)
print(paste0("Now write Year-month: ", yrng[k], " - ", monj, " and use nth of tk: ", nt, " ok"))
}
}
} yk <- read_stars("../data_src/oisst/monthly_anom/202006_anom.nc")
zk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"))
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"))
names(yk)[1] <- "anom"
names(zk)[1] <- "anom"
names(z)[1] <- "anom"
print(range(na.omit(as.data.table(yk)$anom))) #-5.911333 10.598333
print(range(na.omit(as.data.table(zk)$anom))) #-5.39709 16.39870 #old -8.227043 11.727673
print(range(na.omit(as.data.table(z)$anom))) #-7.724047 10.483911 #old -8.361936 6.181935
##gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
##gz <- gplotx(z, "anom", minz = -6, maxz = 16.5, returnx=TRUE)
##gzk<- gplotx(zk,"anom", minz = -8, maxz = 11.0, returnx=TRUE)
#gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
#gzk<- gplotx(zk,"anom", minz = -10, maxz = 7.0, returnx=TRUE)
#gz <- gplotx(z, "anom", minz = -8, maxz = 11.0, returnx=TRUE)
gyk <- plot_anom("../data_src/oisst/monthly_anom/202006_anom.nc", returnx=TRUE)
gzk <- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"), returnx=TRUE)
gz <- plot_anom(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"), returnx=TRUE)layt <- rbind(c(1,1),c(2,2),c(3,3))
grid.arrange(gyk, gzk, gz, layout_matrix=layt)#gyk
#gzk
#gzga <- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2022, "04", "_anom.nc"), #zk,
"anom", #minz = -8, maxz = 11.0,
returnx=TRUE)
gaif (!require("layer")) {
if (!require("remotes")) install.packages("remotes"); library(remotes)
remotes::install_github("marcosci/layer")
library(layer)
}
Reload_sst_nc <- FALSE
if (Reload_sst_nc) {
sst1 <- read_stars(paste0("../data_src/oisst/monthly_sst/202002_sst.nc"))
names(sst1)[1] <- "sst"
sst202002 <- as.data.table(sst1) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
#.[,.(longitude, latitude, heatwave)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst2 <- read_stars(paste0("../data_src/oisst/monthly_sst/201002_sst.nc"))
names(sst2)[1] <- "sst"
sst201002 <- as.data.table(sst2) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst3 <- read_stars(paste0("../data_src/oisst/monthly_sst/200002_sst.nc"))
names(sst3)[1] <- "sst"
sst200002 <- as.data.table(sst3) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst4 <- read_stars(paste0("../data_src/oisst/monthly_sst/199002_sst.nc"))
names(sst4)[1] <- "sst"
sst199002 <- as.data.table(sst4) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
save(sst199002, sst200002, sst201002, sst202002, file="../data_src/stats/sst_1990_2020_stack.RData")
} else {
load("../data_src/stats/sst_1990_2020_stack.RData")
}
tl0 <- layer::tilt_map(sst199002)
tl1 <- layer::tilt_map(sst200002, y_shift=75)
tl2 <- layer::tilt_map(sst201002, y_shift=150)
tl3 <- layer::tilt_map(sst202002, y_shift=225)
map_list <- list(tl0, tl1, tl2, tl3)#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "turbo")zt1 <- fread("../data_src/verify/198201_detrend_anom_yeh.csv")
zt <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 1982, "01", "_anom.nc"))
tt1 <- as.matrix(zt1)
tt <- zt[[1]]
quantile(na.omit(tt1[25,]))## 0% 25% 50% 75% 100%
## -1.57301116 -0.26347324 0.04314169 0.33155985 1.47273052
quantile(na.omit(tt[25,]))## 0% 25% 50% 75% 100%
## -1.866937160 -0.298840664 0.005112544 0.310795255 1.452879071
length(which(is.na(tt1)))## [1] 548586
#[1] 548586
length(which(is.na(tt))) #[1] 474337 !! large diff, should check...## [1] 474337
range(na.omit(as.numeric(tt1)))## [1] -7.042250 3.954143
range(na.omit(as.numeric(tt)))## [1] -7.016757 3.966041
zt <- as.data.table(zt) %>% setnames(3, "anom")
xval <- sort(unique(zt$x))
yval <- rev(sort(unique(zt$y)))
zt1 <- cbind(c(1:1440), zt1)
setnames(zt1, 1, "xid")
zt1a <- melt(zt1, id.vars = "xid", measure.vars = patterns("^V")) %>%
.[,vart:=gsub("V","",variable) %>% as.integer()] %>%
.[,`:=`(x=xval[xid], y=rev(yval)[vart])] %>%
.[,.(x, y, value)] %>% setnames(3, "anom")
setorder(zt, -y, x) ## the same order as original zt data.frame
setorder(zt1a, -y, x)
gz1a <- gplotx(zt1a, "anom", minz = -7.5, maxz = 4, returnx=TRUE)
gzt <- gplotx(zt, "anom", minz = -7.5, maxz = 4, returnx=TRUE)
layt <- rbind(c(1,1),c(2,2))
grid.arrange(gzt, gz1a, layout_matrix=layt)plot(1:720, zt1a[x==121.125,]$anom, type="l", col="blue") #almost the same
points(1:720, zt[x==121.125,]$anom, type="l", col="red", lty=2, add=T)plot(1:720, zt1a[x==10.125,]$anom, type="l", col="green")
points(1:720, zt[x==10.125,]$anom, type="l", col="brown", lty=2, add=T)chk_x <- 360-145+0.125 ## Now try to plot paper 145W, 55N detrend vs anom not detrended
chk_y <- 50.125
plot(1:720, zt1a[x==chk_x,]$anom, type="l", col="black")
points(1:720, zt[x==chk_x,]$anom, type="l", col="magenta", lty=2, add=T)plot(1:1440, zt1a[y==chk_y,]$anom, type="l", col="black")
points(1:1440, zt[y==chk_y,]$anom, type="l", col="magenta", lty=2, add=T)plot(1:1440, zt1a[y==0.125,]$anom, type="l", col="black")
points(1:1440, zt[y==0.125,]$anom, type="l", col="magenta", lty=2, add=T)Reverify_paper <- FALSE
idxt <- which(zt$x==chk_x & zt$y==chk_y) #229821
zt[idxt,] #-0.3095802## x y anom
## 1: 215.125 50.125 -0.3095802
idxm <- idxt %% 1440
idym <- floor(idxt/1440) +1
tt[idxm,idym] == zt[idxt,]$anom ## [1] TRUE
if (Reverify_paper) {
anomx <- c()
anomd <- c()
sstx <- c()
for (i in yrng) {
for (j in 1:12) {
if (i==curryr & j>(currmo-1)) break
monj <- fifelse(j<10, paste0("0",j), paste0(j))
zd <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", i, monj, "_anom.nc"))
anomd <- c(anomd, zd[[1]][idxm, idym])
zk <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
anomx <- c(anomx, zk[[1]][idxm, idym])
zs <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
sstx <- c(sstx, zs[[1]][idxm, idym])
if (i==yrng[1] & j==1) {
testdatex <- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
testdatex<- c(testdatex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
}
save(anomx, anomd, sstx, testdatex, file="../data_src/stats/verify_paper_loc145W50N_sst_anom_1982_202204.RData")
} else {
load("../data_src/stats/verify_paper_loc145W50N_sst_anom_1982_202204.RData")
}
plot(testdatex[1:468], sstx[1:468], type="l", col="magenta")plot(testdatex[1:468], anomx[1:468], type="l", col="magenta", ylim=c(-2,3))
points(testdatex[1:468], anomd[1:468], type="l", col="black", add=T)